#ifndef VECTOR4_HEADER_VIZMF_PROJECT
#define VECTOR4_HEADER_VIZMF_PROJECT

#include <ostream>
#include <cmath>	

#include "vector3.h"

template <typename T>
class Vector4T
{
public:
  typedef T float_type; 
  typedef Vector4T<T> Self;

  Vector4T(float_type a = 0.0, float_type b = 0.0, float_type c = 0.0, float_type d = 1.0)
  {
    x_[0] = a, x_[1] = b; x_[2] = c; x_[3] = d;
  }
  template <typename T1>
  Vector4T(T1 v[4])
  {
    x_[0] = v[0];
    x_[1] = v[1];
    x_[2] = v[2];
    x_[3] = v[3];
  }
  template < typename T2 >
  Vector4T(const Vector3T<T2>& v)
  {
    x_[0] = v[0];
    x_[1] = v[1];
    x_[2] = v[2];
    x_[3] = 1.0;
  }
  Vector3T<float_type> toVector3() const { return Vector3T<float_type>(x_[0]/x_[3], x_[1]/x_[3], x_[2]/x_[3]); }

  float_type length_squared() const { return (x_[0]*x_[0] + x_[1]*x_[1] + x_[2]*x_[2])/x_[3]/x_[3] ; }
  float_type length() const { return std::sqrt( length_squared() ); }
  const Self& normalize() {
    x_[3] *= length();
    return *this;
  }

  float_type& x() { return x_[0]; }
  float_type& y() { return x_[1]; }
  float_type& z() { return x_[2]; }
  float_type& w() { return x_[3]; }
  const float_type& x() const { return x_[0]; }
  const float_type& y() const { return x_[1]; }
  const float_type& z() const { return x_[2]; }
  const float_type& w() const { return x_[3]; }

  const float_type& operator [] (int i) const { return x_[i]; }
  float_type& operator [] (int i) { return x_[i]; }

  const float_type* vertex() const {return x_; }

  template < typename T1 >
  float_type dot_prod(const Vector3T<T1>& v) const { return (x_[0]*v.x_[0] + x_[1]*v.x_[1] + x_[2]*v.x_[2])/x_[3]/v.x_[3]; }
  template < typename T1 >
  Self vec_prod(const Vector4T<T1>& v) const {
    return Self(
      x_[1]*v.x_[2] - v.x_[1]*x_[2], x_[2]*v.x_[0] - v.x_[2]*x_[0], x_[0]*v.x_[1] - v.x_[0]*x_[1],  x_[3]*v.x_[3]
    );
  }
  template < typename T1 >
  bool operator ==(const Vector4T<T1>& other) const
  {
    return x_[0] == other.x_[0] && x_[1] == other.x_[1] && x_[2] == other.x_[2] && x_[3] == other.x_[3];
  }
  template < typename T1 >
  Self operator + (const Vector4T<T1>& other) const
  {
    return Self(x_[0]/x_[3] + other.x_[0]/other.x_[3], x_[1]/x_[3] + other.x_[1]/other.x_[3], x_[2]/x_[3] + other.x_[2]/other.x_[3], 1.0);
  }
  Self operator - () const
  {
    return Self(-x_[0], -x_[1], -x_[2], x_[3]);
  }
  template < typename T1 >
  Self operator - (const Vector4T<T1>& other) const
  {
    return *this + (-other);
  }
  template < typename T1 >
  Self operator += (const Vector4T<T1>& other)
  {
    return *this = *this + other;
  }
  template < typename T1 >
  Self operator -= (const Vector4T<T1>& other)
  {
    return *this = *this - other;
  }
  Self operator * (const float_type& a) const
  {
    return Self( a * x_[0], a * x_[1], a * x_[2], x_[3]);
  }
  Self operator / (const float_type& a) const
  {
    return Self( x_[0] / a, x_[1] / a, x_[2] / a, x_[3]);
  }
  Self operator *= (const float_type& other)
  {
    return *this = *this * other;
  }
  Self operator /= (const float_type& other)
  {
    return *this = *this / other;
  }
  template < typename T1 >
  Self perpendicular(const Vector4T<T1>& other) const
  {
    return this->vec_prod(other);
  }
  Self perpendicular() const
  {
    Self result = perpendicular(Self(0.0, 0.0, 1.0, 1.0));
    if (result.length_squared() < ZERO) return result; else return perpendicular(Self(1.0, 0.0, 0.0, 1.0));
  }
  static Self X()
  {
    return Self(1.0, 0.0, 0.0, 1.0);
  }
  static Self Y()
  {
    return Self(0.0, 1.0, 0.0, 1.0);
  } 
  static Self Z()
  {
    return Self(0.0, 0.0, 1.0, 1.0);
  }
  static Self W()
  {
    return Self(0.0, 0.0, 0.0, 1.0);
  }
  template <typename T1, typename T2>
  static int compare(const Vector4T<T1>& lhs, const Vector4T<T2>& rhs)
  // returns -1 if lhs < rhs, 0 if lhs == rhs and 1 if lhs > rhs
  {
    if (fabs(lhs.x_[0] - rhs.x_[0]) < ZERO)
    {
      if (fabs(lhs.x_[1] - rhs.x_[1]) < ZERO)
      {
        if (fabs(lhs.x_[2] - rhs.x_[2]) < ZERO)
        {
          return 0;
        } else {
          if ((lhs.x_[2] - rhs.x_[2]) < 0.0) return -1; else return 1;
        }
      } else {
        if ((lhs.x_[1] - rhs.x_[1]) < 0.0) return -1; else return 1;
      }
    } else {
      if ((lhs.x_[0] - rhs.x_[0]) < 0.0) return -1; else return 1;
    }
  }

private:
  float_type x_[4];
};

template < typename T >
inline Vector4T<T> operator * (const double& a, const Vector4T<T>& vec)
{
  return vec * a;
}

#ifdef DEBUG

template < typename T >
inline std::ostream& operator << (std::ostream& strm, const Vector4T<T>& vec)
{
  return strm << "(" << vec[0] << ", " << vec[1] << ", " << vec[2] << ", " << vec[3] << ")";
}

#include <QDebug>

template < typename T >
inline QDebug& operator << (QDebug& strm, const Vector4T<T>& vec)
{
  return strm << "(" << vec[0] << ", " << vec[1] << ", " << vec[2] << ", " << vec[3] << ")";
}

#endif // DEBUG

typedef Vector4T<float> Vector4f;
typedef Vector4T<double> Vector4d;
typedef Vector4f Vector4;

#endif // VECTOR4_HEADER_VIZMF_PROJECT
